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Abstract 

Structural correlations at a liquid-solid interface were explored with molecular dy- 
namics simulations of a model aluminium system using the Ercolessi-Adams po- 
tential and up to 4320 atoms. Substrate atoms were pinned to their equilibrium 
crystalline positions while liquid atoms were free to move. The density profile at the 
interface was investigated for different substrate crystallographic orientations and 
temperatures. An exponential decay of the density profile was observed, p(z) ~ e~ KZ , 
leading to the definition of k as a quantitative measure of the ordering at the liquid 
solid interface. A direct correlation between the amount of ordering in the liquid 
phase and the underlying substrate orientation was found. 
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1 Introduction 



Metal-ceramic interfaces play a prominent role in a variety of technological 
applications and processes that range from electronic devices, protective coat- 
ings, high-temperature structural components, and liquid phase joining pro- 
cesses. The functionality of these systems depends crucially on their macro- 
scopic properties such as fracture strength, yield, and electrical conductivity. 
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These properties are strongly correlated with microscopic details of the metal- 
ceramic interface such as wetting, chemistry, diffusion, and structure. Corre- 
lating macroscopic properties to the structure and chemistry of interfaces is 
one of the most intriguing topics in materials science. Experimental studies 
of the atomistic structure of a solid- liquid internal interface are technically 
difficult to conduct [1,2]. Therefore, atomistic simulations of metal ceramic 
interfaces can serve as an important tool to understand and predict the effect 
of the interface region on the material properties. Ab-initio electronic structure 
calculations provide detailed information of the chemical bonding across the 
interface, but are very expensive in terms of computer power and time. Conse- 
quently, they are limited to modeling small systems, which are not of sufficient 
size to contain large structural defects. Moreover, a dynamical simulation of 
a liquid-solid interface is beyond both system size and CPU limitations for 
ab-initio calculations. 

Atomistic simulations, such as Molecular Dynamics (MD) or Monte Carlo 
permit the controlled study of these systems at the atomistic level for a large 
number of atoms and for large structures. However the main limitation to such 
simulations is the lack of appropriate interatomic potential schemes which 
can model both metallic and ionic bonding across the interface. Nevertheless, 
simplified models can be used to obtain qualitative basic insights into the 
problem. In this study we introduce a model system in which the ceramic is 
assumed to be composed of atoms pinned to their equilibrium lattice positions, 
while the metal atoms are free to evolve under the influence of their interatomic 
potential. 

The atoms of a liquid metal which are adjacent to a rigid crystalline substrate 
are in an environment which is strongly affected by the symmetry of the 
underlying substrate. Theoretical studies [3], which are mainly computational, 
have shown that ordering occurs in the first layers of the liquid adjacent to 
the crystal surface. The same result emerges from experimental studies of 
solid-liquid interfaces [1,3]. A particularly interesting study with respect to 
the current project was conducted by Huisman et al.[l], who investigated the 
interface structure of liquid gallium in contact with a diamond (111) surface. 
They observed pronounced layering of the liquid metal density profile which 
decays exponentially with increasing distance from the diamond. Moreover, 
the interlayer spacing was equal to the repeat distance of (100) planes of 
upright gallium dimers in solid a-Ga. Thus it appears that the liquid near the 
surface assumes a solid-like structure similar to the a-phase. 

Consequently, at a solid-liquid interface, there should be a transition from a 
solid-like phase in the liquid near the interface to a liquid like phase further 
away from the interface. The extent over which such a transition occurs should 
be directly related to the structure of the underlying surface. In the current 
research, simulations have been conducted to study the density profile and 
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structure of the liquid-metal/hard-wall interface as a function of temperature 
and substrate structure. In particular we show that the decay of the density 
profile is quantitatively and qualitatively related to the underlying structure 
of the substrate. 



2 System and Simulation Method 

The initial configuration for the simulations was a crystalline face-centered- 
cubic (fee) aluminium lattice. A liquid-solid interface was formed by pinning 
the atoms in the first few layers to their ideal lattice positions, and allowing the 
rest of the atoms to move freely under the effect of the interatomic potential. 
The lattice parameter of the fee aluminium crystal was a = 4.1 A, which is 
close to the value obtained for aluminium at the melting point (T = 940K) 
with the potential model used here [4] (see below for more details of the 
potential). Three orientations for the terminating substrate plane were used: 
(100), (110), and (111). An example of a simulation cell with a (110) interface is 
shown in Fig. 1. The light gray atoms were fixed to their equilibrium positions 
throughout the simulation. The number of particles in each sample, as well as 
the geometry, are described in Table 1. 

In the plane of the interface (the xy plane) periodic boundary conditions 
were applied. In the direction perpendicular to the interface (the z direction), 
the boundary conditions are expected to simulate the bulk media on either 
side of the interface. On the rigid (substrate) side where atoms are fixed to 
crystalline positions, it is sufficient to require that the extent of the region in 
the z direction to be larger than the cutoff of the interaction potential. In this 
way liquid atoms near the interface do not "see" the bottom of the rigid layer, 
hence it acts like a semi infinite bulk system. On the liquid (metal) side the 
situation is more difficult. If periodic boundary conditions are applied, the first 
rigid layer (see Fig. 1) and the last liquid layer will interact forming a second 
interface in the system. These two interfaces can then interact with each other, 
unless the dimension in the z direction is very large so that the interaction 
is negligible. In addition, due to the confinement of the simulation cell, the 
system will not be able to respond to volume changes caused by stress at the 
interface. One way to overcome this difficulty is to use a model in which the 
sample is enclosed by two moving rigid walls as proposed by Lutsko et al. [5]. 
Here, to allow for volume deformations, the z dimension is considered as a 
dynamic variable and a formulation is introduced to allow a flexible length of 
the sample normal to the interface in response to stresses in the simulation 
cell. 

In this study a liquid layer is deposited above the solid layer, and then a 
vacuum region is inserted with periodic boundary conditions in all directions. 
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In this way we have one free liquid interface and one internal solid-liquid 
interface. Provided that the height of the liquid layer is large enough in the z 
direction, there will be no interaction between the free liquid surface and the 
internal rigid-liquid interface. In addition the system will be free to respond 
to stresses at the interface since there is nothing to limit the liquid from 
expanding in the z direction. 

Since periodic boundary conditions were also applied in the z direction, there 
are two requirements on the extent of the vacuum region. Firstly, it should 
be larger than the cutoff distance of the interatomic potential, otherwise the 
upper most liquid layer will interact with the bottom of the rigid layer. Sec- 
ondly, if liquid atoms evaporate from the free surface and enter into the vac- 
uum region, there is a chance that these vapor atoms would reach the end of 
the vacuum and interact with the bottom region via the periodic boundary 
condition. However, as discussed below, there are practically no vapor atoms 
during the simulation time. Hence it is sufficient to have a large enough vac- 
uum based on the cutoff distance only. In practice we chose a vacuum region 
of 15 times the extent of the cutoff distance. The system was simulated us- 
ing an MD technique, which consists of the numerical integration of Newton's 
equation of motion for the various atoms [6]. The velocity- Verlet integration 
algorithm [6,7] was used in the simulations. For the interatomic potential, an 
embedded atom potential developed by Ercolessi and Adams (EA) [8] was 
used. This potential was constructed by the so called force matching method, 
whereby the potential was fitted to a very large amount of data obtained from 
both experiment and first principle calculations, with emphasis given to match 
the interatomic forces obtained from the potential to those obtained from first 
principles. The potential has been tested in detail for aluminium[4,9], and was 
found to be consistent with experimental results. As an example, the calcu- 
lated melting point for aluminium is T=939±5K in excellent agreement with 
the experimental value of T=933.6K. 

One advantage for using a more realistic metallic potential (instead of a simpler 
pair potential such as a Leonard- Jones (LJ) 6-12 potential), is that it gives the 
realistic low vapor pressure characteristic of liquid metals. Simulations that 
were performed with a LJ potential resulted in evaporation of a large portion 
of the liquid. In the case of liquid aluminium, the low vapor pressure results in 
as little as 10~ 8 vapor atoms in a simulation box of cubic side equal to ~ 40 A, 
therefore, in effect no vapor atoms are observed during the simulations. 

As a validation of our implementation of the EA potential we calculated the 
calorie curve for a bulk liquid system i.e. the internal energy as a function 
of temperature E(T), the pressure, and the translational order parameter as 
functions of temperature. In all cases we observed a discontinuous jump at a 
temperature of T = 930±15K, which is consistent with the calculated melting 
temperature of Ercolessi and Adams[8]. Of course this is not the best way to 
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measure the melting point, but it suffices as a preliminary check of our correct 
implementation of the EA potential in our code scheme. 



As stated above, the size of the substrate (the rigid layers) is determined 
by the cutoff of the interaction potential. The range of effective interactions 
in an EAM type of potential is actually twice the cutoff radius of the glue 



interact with the bottom of the rigid sample, which is effectively a semi-infinite 
bulk. A simple way of modifying the substrate liquid interaction is to reduce 
the number of substrate planes from this value, and we have included such 
results here. 

The atoms in the fixed region were excluded from the equations of motion, 
although the forces they exert on the adjacent layers were included. In this 
way these fixed layers can be thought of as being a part of a different material 
with a much higher melting temperature, such as a ceramic. 

The temperature of the simulations was controlled by a simple ad-hoc rescaling 
of the velocities of the particles, so that the required average kinetic tempera- 
ture was reached [7]. Two schemes of temperature control were used; rescaling 
of all particle velocities, and rescaling only of the velocities of atoms adjacent 
to the fixed region, i.e. of the two layers next to the fixed one. The time step 
of the MD integration in normalized units was t = 0.02, where the normalized 
unit is t = 4.25 x 10~ 14 in real seconds. 



2.1 Computing the Density Profile 

The density profile p(z) is defined as the average density of particles in a slice 
of width Az parallel to the hard wall surface and centered around x. The 
simulation cell is divided into equal layers or bins parallel to the interface. 



1 Normally the d-spacing of a family of planes with Miller indices (hkl) is given by: 
dhki = aiatt/Vh 2 + k 2 + I 2 - For example in the case of (110) planes: duo = ^latt/V^- 
However, this is not the smallest distance between two consecutive planes having 
the same in-plane structure as the (110) plane. The family of planes that includes 
all planes with a (110) structure is in fact the (220) family. Hence in calculating 
the d-spacing for the (110) surface, we used cfeo- Similarly for the (100) family, we 
used d 2 oo- 
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The expression for the density profile is 

(N z ) 



Piz) = 



L x L y Az 



where L x and L y are the x and y dimensions of the cell, respectively, and z 
is perpendicular to the interface, Az is the bin width, and N z is the number 
of particles between z — Az/2 and z + Az/2 at time t. The angled brackets 
indicate a time average. 

In order to reduce the statistical error of the sampling, a proper choice of bin 
width must be made. Very small bin widths results in too few particles at each 
time step, hence a large scatter of the data. Very large bins will not show the 
actual dependence of the density profile over distance. Two basic width scales 
have been used: a coarse scale, in which the width of the bins was set equal 
to the bulk crystal d-spacing for a particular orientation, and a fine scale for 
which each coarse scale bin was divided into 10 or 25 sections. 



2.2 Equilibration and Computation of Averages 



For each system described in Table 1 a series of simulations was conducted at 
temperatures ranging from T = 800K up to T = 2000K in steps of 50K. To 
obtain an equilibrated sample at a particular temperature, an ideal configura- 
tion is annealed by heating the system in steps of 50K to a higher temperature 
than required, and then cooling down in steps of 50K to the target tempera- 
ture. Fig. 2 shows one such annealing schedule for system No. 1 (see Table 1). 
After annealing, the system is further equilibrated at the target temperature 
for at least 50, 000 MD steps for the smaller systems and 100, 000 steps for 
the larger systems, before sampling begins. 

Measurements of the density profile are obtained by accumulating p{z) via 
equation (1) over 20,000 to 50,000 MD steps, then averaging over these to 
produce a single measurement. In equilibrium, these measurements of p{z) 
do not change substantially in time, and are independent of the annealing 
schedule. 



3 Results and Discussion 



The density profiles p(z) of the liquid part of the samples are shown at two 
temperatures in Figs. 3, 4, and 5 for the (111), (HO), and (100) interfaces, 
respectively. The inserts in the figures are enlargements of the profiles of the 
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first few layers, with the abscissa given as the average number of particles per 
bin (N z ) [see equation (1)], rather than as the density. These density profiles 
show large oscillations corresponding to the layering of the liquid near the 
hard wall. The oscillations dampen gradually within the interfacial region, 
until the uniform density of the liquid phase is reached, which is the same in 
all systems: p\ = 0.0051 ± 0.0005 atoms/A 3 . As expected, the density profile 
at a higher temperature decays faster as a function of the distance from the 
interface than the density profile at a low temperature. 

In order to facilitate a direct comparison between the density profiles, the den- 
sity is normalized by the magnitude of the first peak as shown in Fig. 6, where 
the density profiles for these systems are plotted at T = 1000K. It is easy to 
see from Fig. 6 that ordering at the (111) and (100) interfaces extends further 
into the liquid than at the (110) interface. However, the number of density 
peaks is very similar in all directions, while the distance between consecutive 
peaks varies substantially with interface direction. This distance, or interlayer 
spacing can be identified with the d-spacing of the quasi solid region within 
the liquid. The calculation of the d-spacing is made with an algorithm that 
searches for maxima in the density profile, and then calculates the difference in 
the locations of each consecutive maxima, giving the d-spacing. The d-spacing 
is plotted in Fig. 7 at T = 1000K. The data in Fig. 7 includes the analysis 
of the d-spacing even well inside the liquid region, where the "maxima" are 
in random locations and are due to noise in the data. This explains the large 
scatter of data at a distance from the interface, and can serve as an indication 
of the point at which the ordering terminates. For example, it is observed that 
the ordering of the (111) interface extends up to almost 25 A, before becoming 
too noisy. It is also clear that the interlayer spacing in the quasi-liquid is the 
same as that in the rigid substrate, at least in the first few layers. However, 
the (100) layers show a very large expansion of the d-spacing in the direction 
of that of the (111) system. Similar behavior was previously observed [10-12], 
and arises from the fact that perturbations within the liquid are more energet- 
ically favored if their periodicity is 27i/Qq, where Qo is the wave vector at the 
first peak of the structure factor S(Q). For an fee liquid metal this distance is 
very close to dm [12]. The d-spacing starts to deviate from the bulk value far 
away from the substrate, as can be seen for the (111) case, where d alternates 
between two values very close to rf m before becoming completely random. For 
the (110) case a shift in the interlayer spacing is not observed. This may be 
attributed to the large difference between duo and dm. 

In order to analyze the decay of the density profile as a function of temperature 
and substrate orientation, the envelope of the density profile was extracted and 
plotted against distance from the surface, as shown in Fig. 8 (a). The profile is 
then fitted assuming an exponential decay: p{z) ~ exp (—kz), where k — l/£, 
and £ is the correlation length at the interface. The parameter k can then be 
used to quantitatively describe the amount of disorder at the interface. The 
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exact function that has been used in the fitting has an extra constant term, b, 
to account for the background liquid density, and a normalization factor, a: 

p( z ) = ae~ KZ + b (2) 



This form of decay of the density profile is typically obtained from a mean 
field treatment of binary fluid interfaces [12,13]. A good check of the validity 
of this assumption is that log(p(z) — b) is linear in z. A few typical examples 
are shown in Fig. 8 (b). The large scatter for points far from the solid-liquid 
interface, i.e. well within the normal liquid phase, is due to the sensitivity of 
the logarithm to small arguments, since the difference between p(z) and b for 
points inside the liquid is very small. 

The results of the density profile analysis for the three interfaces, namely 
(111), (110) and (100) are shown in Fig. 9. Here the disorder parameter k is 
plotted against temperature for the three orientations (111), (HO) and (100), 
for various systems as indicated in the figure. It is clear that the finite size 
effects are very small, as can be seen from the compatible k values for systems 
No. 1 and 5 (see Table 1) and systems No. 2 and 3. Moreover, we note that the 
width of the underlying rigid substrate does not alter the qualitative behavior 
of k, as seen from the results for systems No. 6 and 7. However, there is 
a tendency for k to be larger for a thick substrate, as evident mainly from 
the results for system 7. This can be attributed to the "gluing" nature of 
the potential, namely that the energy of an atom is lower when it has more 
neighbors. Hence when there are more substrate layers the energy of the overall 
system is lowered, and the net force exerted on the liquid is lower than with 
a thinner substrate which leads to less ordering in the liquid phase. 

The amount of disorder, k, is consistently smaller for the (111) and (100) 
surfaces than for the (110) system, as shown in Fig. 9. In other words, the 
ordering at the (111) and (100) interfaces is larger, and extends further into 
the liquid, than the ordering at the (110) interface. 

The close-packed surfaces (100) and (111) in an ideal system have a lower 
energy than the (110) surface, and therefore are more stable against thermal 
effects. The formation energies of an ideal surface, as calculated with the EA 
potential [9], and with LDA [14] are shown in table 2. The formation energies 
of an ideal (111) and (100) surfaces are similar, with the difference between 
them being only ~ 0.1 eV/atom with the EA potential. The formation energy 
of the (110) surface is 0.26 eV/atom higher than that of the (100) surface. 
The similarity of the formation energies of the ideal (111) and (100) systems 
is compatible with the behavior of k in Fig. 9, where k for the (111) and 
(100) systems are almost indistinguishable. This suggest that the the energy 
balance between the different surfaces is not altered by the presence of a liquid 
phase, although this is not immediately obvious, and explicit calculations have 
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to be made for the different energies of a solid/liquid interface. However, the 
behavior of k strongly suggests that the energy balance is not drastically 
changed. 

Fig. 9 also shows that the rate of increase of disorder is larger in the [110] 
direction than for the closed packed surfaces (100) and (111). Again this is 
consistent with the lower formation energy of ordered (111) and (100) surfaces. 
Another striking difference is that as temperature increases the amount of dis- 
order in the closed packed directions seems to saturate while that for the (110) 
surface seems to increase approximately linearly with temperature. There is 
even a slight decrease in k i.e. an increase of ordering, at high temperatures. 
Note that data for high temperatures (T > 1400/T) should be analyzed with 
caution since at these temperatures there are one or two peaks at most and 
the fitting procedure is no longer accurate. 

To summarize, the results indicate that the interaction between the ordered 
solid and the liquid induced layering oscillations within the liquid with a pe- 
riodicity of the interlayer spacing d. If this periodicity is also commensurate 
with the natural periodicity of the liquid, as given by 2n/Q (see above), then 
these oscillations are enhanced, and survive a further distance into the liquid. 
This is the case for the (111) interface. In the case where the d-spacing is 
very close to 2n/Q , such as for the (100) interface, then after some distance 
(about three layers in Fig. 7) when the effect of the hard-wall has faded, the 
d-spacing switches to a value close to 2ir/Q . When the d-spacing is far from 
2tt/Qo, then the oscillations fade away after a short distance. 

This suggests that k depends strongly on the d-spacing of the underlying sub- 
strate. In fact a reasonable order of magnitude estimate for Kjfi clS db function of 
orientation n, can be obtained by assuming that to first order, the orientation 
dependence is related only to the d-spacing by: 

K ~ \ ^ 

With the value used for the lattice parameter we obtain: 

Klu = 0.433A" 1 , «no = 0.689A" 1 , and k 100 = 0.488A" 1 . (4) 

Inspection of Fig. 9 shows that this occurs for an intermediate temperature 
of approximately 1200K. Again we obtain similar k values for the (111) and 
(100) surfaces, even when neglecting the jump in d 1Q0 which would make k 
even more comparable. 

The assumption that the layering periodicity within the liquid follows that 
of the underlying substrate allows us to explain the different rates at which 
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disorder increases with temperature for the different orientations, as shown in 
Fig. 9. If the d-spacing in the liquid region adjacent to the substrate is large, 
as the case for the closed-packed surfaces (111) and (100), then the overlap 
between the density peaks is very small, as can be seen from Fig. 6. For open 
surfaces such as (110) the d-spacing is smaller and the overlap is larger. A small 
overlap leads to a small number of atoms in between the peaks, as is clearly 
observed from the deep valleys in Fig. 6, while a large overlap leads to a higher 
density of atoms in between the peaks. In a metallic system, the energy of an 
atom is lower when it is in a dense environment, and as a result atoms tend 
to reside in a dense layer rather than in a less dense valley. Breaking the order 
of the layers essentially means that more atoms leave the ordered layers and 
reside between the layers. Therefore, in orientations with a large d-spacing, 
where valleys are almost void of atoms, the barrier to such interlayer diffusion 
is large. So for close-packed surfaces (111) and (100), where order prevails 
up to higher temperatures, the peaks remain well separated, and hence as 
temperature increases, there is almost no increase in k. For the (110) interface, 
the larger overlap between the peaks allows for easy interlayer diffusion, and 
hence k increases readily. 

We presume that the link between the d-spacing of the liquid and the substrate 
is not causal, but is due to the imposition of a parallel structure on the liquid 
by the substrate surface. If the liquid atoms under the first peak in the liquid 
density are strongly correlated with the surface atomic structure, then the fact 
that our substrate was chosen to have nearly the same atomic density as the 
liquid ensures that the distance to the second liquid peak is the same as the 
interplanar spacing in the substrate. It might be interesting to demonstrate 
this by varying the d spacing of the substrate without changing the density 
within the lattice planes parallel to the surface. We would not expect this to 
lead to significant changes in the liquid peak separations. Other extensions of 
this study would be to introduce a substrate with a different lattice structure 
to that of the frozen liquid, eg. bcc in the case of liquid aluminium. Finally, 
the effect of varying the lattice parameter of the substrate, such that it is no 
longer commensurate with the frozen liquid in its normal crystal structure, 
may induce some interesting effects for further study. 



4 Summary and Conclusions 

Atomistic simulations on a model liquid-solid interface have been performed, 
and the density profile of the liquid phase as a function of temperature and 
substrate orientation has been studied. The substrate was composed of alu- 
minium atoms pinned to an fee lattice structure, while the liquid metal was 
composed of normal liquid aluminum. Substantial ordering was found in the 
liquid near the interface. The density profiles showed large oscillations cor- 



10 



responding to the layering of the liquid near the interface. The oscillations 
dampen gradually within the interfacial region according to p(z) ~ e~ KZ , until 
the density of the liquid phase is reached. The ordering was found to ex- 
tend further into the liquid region at close-packed interfaces such as (111) and 
(100). Moreover, the periodicity of the liquid layering was found to be strongly 
correlated to the inter-planar spacing of the rigid substrate, which in turn is 
related to the planar structure of the underlying substrate. 

The decay parameter k was defined as a disorder parameter which gives the 
extent at which the layering extends into the liquid. It was found that k is con- 
sistently larger for the (110) interface, than for the (111) and (100) interfaces. 
Moreover, as temperature increases, k for the open-surface (110) increases ap- 
proximately linearly with temperature, while for the (111) and (100) interfaces 
it saturates. The interlayer spacing in the liquid also determined the amount 
of overlap between the density peaks. The smaller the d-spacing, the larger 
the overlap between the ordered liquid layers, which leads to an enhanced 
interlayer diffusion resulting in a faster increase of disorder. 
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Figure 1: An example of a cell for the simulation of a system with a free (110) surface. 

The total number of atoms is 1944, out of which 432 atoms are 'solid' (light 
gray) and the rest are "liquid" (dark gray) aluminium atoms. 

Figure 2: An example of an annealing schedule for system number 5 (see Table 1). 

The state of the system at the final (target) temperature is compared with 
the initial state. The schedule is repeated until the system is in equilibrium, 
that is when the energy and density profile in the final state and initial state 
are the same. 

Figure 3: The density profile p(z) as a function of distance from the hard wall z, 
is shown for a (111) system at two temperatures: 7\ = lOOOi^, and T 2 = 
12007^. The insert is a plot of the average number of atoms per bin (N x ) 
for the region near to the rigid wall. Data is from a simulation on system 
number 4 described in Table 1. 

Figure 4: The density profile p(z) as a function of distance from the hard wall z, 
is shown for a (110) system at two temperatures: 7\ = lOOOif, and T 2 = 
1200-fT. The insert is a plot of the average number of atoms per bin (N x ) 
for the region near to the rigid wall. Data is from a simulation on system 
number 2 described in Table 1. 

Figure 5: The density profile p{z) as a function of distance from the hard wall z, 
is shown for a (100) system at two temperatures: Ti = 1000/T, and T 2 = 
1200i^. The insert is a plot of the average number of atoms per bin (N x ) 
for the region near to the rigid wall. Data is from a simulation on system 
number 1 described in Table 1. 

Figure 6: Normalized fine scale density profile for (a): (111), (b):(110), and (c) (100) 
interfaces at the same temperature (1000K). 

Figure 7: Layer separation across the samples at T=1000K, the dark symbols are 
within the rigid substrate, and the open symbols are within the liquid. 

Figure 8: The envelope of the density profile obtained for a (110) sample at T = 
lOOOi^. System numbers are defined in Table 1. 

Figure 9: The disorder parameter k vs. temperature for the three surfaces. 
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Table 1 



Simulation cell geometry. 



System No. 


Interface plane 


No. of planes in each direction 


No. of rigid planes 


Number of atoms 


1 


(100) 


22 [ioo] x 14[oio] x 14[ooi] 


3 


2156 


2 


(110) 


12[m] x 30[no] x 36[2n] 


6 


2160 


3 


(110) 


18[ m ] x 30[ 110 ] x 48[ 2 n] 


6 


4320 


4 


(111) 


18[ m ] x 18[ 110 ] x 36[ 2 n] 


3 


1944 


5 


(100) 


40[ioo] x 14[oio] x 14 [ooi] 


3 


3920 


6 


(111) 


32[ U1 ] x 18[ 110 ] x 36[ 2 n] 


14 


3456 


7 


(110) 


12[ U1 ] x 34[ 110 ] x 36[2n] 


10 


2448 



Table 2 

Ideal surface formation energies for aluminium as calculated with the Ercolessi- 
Adams potential 

System 



9] and with LDA [14]. 
EA [eV/A 2 ] [eV/atom] 



(111) 
(100) 
(110) 



0.054 
0.059 
0.065 



0.38 
0.48 
0.74 



LDA [eV/A 2 ] [eV/atom]] 



0.07 
0.071 
0.08 



0.48 
0.56 
0.89 



[110] 




Fig. 1. An example of a cell for the simulation of a system with a free (110) surface. 
The total number of atoms is 1944, out of which 432 atoms are 'solid' (light gray) 
and the rest are "liquid" (dark gray) aluminium atoms. 



14 



1700 



1600 



o 1500 

Q. 



1400 - 



1300 - 



1200 

500 1000 1500 2000 

MD steps x 10 3 



Fig. 2. An example of an annealing schedule for system number 5 (see Table 1). The 
state of the system at the final (target) temperature is compared with the initial 
state. The schedule is repeated until the system is in equilibrium, that is when the 
energy and density profile in the final state and initial state are the same. 
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Position along the [1 1 1] direction x, in [A] 

Fig. 3. The density profile p(z) as a function of distance from the hard wall z, is 
shown for a (111) system at two temperatures: T\ = 1000-fT, and T2 = 1200-fT. The 
insert is a plot of the average number of atoms per bin (N x ) for the region near to 
the rigid wall. Data is from a simulation on system number 4 described in Table 1 
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Position along the (110) direction x, in [A] 



Fig. 4. The density profile p(z) as a function of distance from the hard wall z, is 
shown for a (110) system at two temperatures: T x = 1000K, and T 2 = 1200^. The 
insert is a plot of the average number of atoms per bin (N x ) for the region near to 
the rigid wall. Data is from a simulation on system number 2 described in Table 1. 
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Position along the [100] direction x, in [A] 

Fig. 5. The density profile p{z) as a function of distance from the hard wall z, is 
shown for a (100) system at two temperatures: T\ = 1000JT, and T<i = 1200-fT . The 
insert is a plot of the average number of atoms per bin (N x ) for the region near to 
the rigid wall. Data is from a simulation on system number 1 described in Table 1. 
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Fig. 6. Normalized fine scale density profile for (a): (111), (b):(110), and (c) (100) 
interfaces at the same temperature (1000K). 
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Fig. 7. Layer separation across the samples at T=1000K, the dark symbols are 
within the rigid substrate, and the open symbols are within the liquid. 
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Fig. 8. The envelope of the density profile obtained for a (110) sample at T = 1000-ftT. 
System numbers are defined in Table 1. 
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9. The disorder parameter k vs. temperature for the three surfaces. 



22 



